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Abstract 

We present an algorithm for generating Poisson-disk 
patterns taking 0{N) time to generate N points. The 
method is based on a grid of regions which can con- 
tain no more than one point in the final pattern, and 
uses an explicit model of point arrival times under a 
uniform Poisson process. 



1 Introduction 

There is a long-standing interest in Poisson-disk 
patterns in the graphics community, primarily for 
their use in sampling |Yel83| ICooSBj IMit87j . There 
have been many algorithms for generating such pat- 
terns. Direct implementation of "dart-throwing" 
|Mit87[ IMF92) produces true Poisson-disk patterns, 
but is slow to converge. Approximations from re- 
laxation |Llo83| or tiling pDJ04, HDKQl] can pro- 
duce patterns similar to Poisson-disk patterns more 
efficiently. Recently, exact methods taking log-linear 
time {0{N log N) where N is the total number of 
points) have been described |DH06[ IJon06| . as well 
as a method with empirical 0{N) speed, but lacking 
a rigorous proof of this performance WCE07. 

We present an algorithm with provable 0{N) per- 
formance. The algorithm maintains two data struc- 
tures: a grid of regions in which points might still 
be inserted, and a bucket (i.e., an unordered set) of 
regions where a point will be generated (a subset of 
the grid). At each step of the algorithm, a region 
is taken from the bucket, a new point is inserted in 
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that region, and nearby regions are updated and pos- 
sibly added to the bucket. The bucket is only empty 
when no more points can be added (i.e., the pattern 
is maximal). The work for each iteration is 0(1), for 
a total cost of 0{N). 

Python source code is included in the ancillary data 
with this paper. 



2 Background 

The Poisson Disk distribution can be defined as the 
limit of a uniform two-dimensional Poisson process 
with a minimum-distance rejection criterion. Succes- 
sive points are independently drawn from the uniform 
distribution on [0, 1]^. If a new point is at least dis- 
tance R from all points already accepted, it is also 
accepted. Otherwise, it is rejected. We call this the 
naive algorithm. The choice of R controls the min- 
imum distance between points (for N points in the 
unit square, ttR^N/A « 0.548 as i? -> |DWJ91| V 

Efficient algorithms for Poisson-disk patterns rely 
on generating new points in regions where they 
are guaranteed (or highly probable) to be accepted 
|DH06l iJoKoel IWCEOT] . In order to guarantee equiv- 
alence of results with the naive algorithm, these 
methods have used O(logA^) area- weighted binary 
search to find where to insert a new point |DH061 
IJon06) . or weighted spatial indexing jWCE07| with 
theoretical O(logA^) but empirical 0(1) cost. 



3 Method 

Our algorithm can be seen as an optimization of the 
naive algorithm using a spatial data structure. We 




Figure 1: A point p is shown with its neighbor grid 
squares in grey. If p's arrival time is earher than any 
of its neighbors, it wih be accepted and added to the 
output. The free regions of p's neighbors wiU then 
be updated to the dark gray areas. This may result 
in points such as q being invalidated. In g's case, a 
replacement point q' will be generated in its new free 
region, Aq> . The time of arrival of q' will be ig + t+ , 
where i_|_ is drawn from an exponential distribution 
parameterized by the size of the updated free region. 



store a grid with spacing < i?/\/2 such that no more 
than one point can land in any grid square in the final 
pattern. We model (implicitly) a uniform 2D Poisson 
process on [0, 1]^, with rate A = 1, by storing at each 
grid square the location and time of the earliest point 
landing in that square. These points and their arrival 
times will be updated as the algorithm progresses. 

Each grid square has three associated pieces of 
data: the free region within that square where new 
points might be generated, a random point within 
that free region, and a time-of-arrival for that point 
under a the Poisson process. Initially, the free region 
for each grid square is the entire square, each point 



is chosen uniformly within its square, and the times 
of arrival are drawn from Aoe~ °*, where Aq is the 
area of a grid square. 

We also define a neighbor relationship from points 
to grid squares, where the neighbors of a point are 
any grid squares within R of the point (see figure [I]) . 

The first insight of our paper is that any point p 
that has time-of-arrival tp lower than any of its neigh- 
bors can be added to the output immediately, as this 
indicates that p arrives before any other point that 
could prevent it from being accepted. On acceptance, 
the free regions of p's neighbors are updated (see fig- 
ure [l]). 

It is possible that accepting p will invalidate a point 
q from another grid square with ||p — 9II2 < R and 
tq > tp, in which case q is removed from the grid and 
a new point g' in the updated free region is created 
with a new and later tqi (see figure [l]) . 

The second key insight of our algorithm is that the 
new tgi should be tg plus a random variable drawn 
from the exponential distribution parameterized by 
the area of the updated free region Aq>, i.e., tq> = 
tq + t+, where t+ is drawn from Aq>e~'^i''^+ . 

The logic is as follows. The points p and q represent 
the first arrivals in their respective grid squares, at 
times tp and tq (in the naive algorithm). When q is 
invalidated by p, the time until another point arrives 
(in the now smaller free region) is modeled by an 
exponential process with parameter Aq> . 

To track which points are candidates for accep- 
tance, we traverse the grid and identify every point 
that has a time of arrival earlier than any of its neigh- 
bors (ignoring neighbors that have already had their 
points accepted). We term these points locally early, 
and add them to a bucket (an unordered set). At each 
iteration, we can take any point from the bucket, and 
add it to the output pattern. 

Accepting p may lead to new points becoming lo- 
cally early, which are then added to the bucket. Like- 
wise, if a point q is invalidated by p's acceptance, 
points with g's grid square in their neighbors may 
become locally early, as g's replacement q' will have 

Iq' -> Tq. 

Each iteration is 0(1) provided we can update, 
compute the area of, and sample uniformly from the 
free space of a grid square in 0(1) time. Previous 
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Figure 2: Performance (PD samples per second) over 
a range of exclusion radii (0.64 to 0.004 in a geo- 
metric progression with ratio 0.8) within the unit 
square. The number of samples generated by the 
Poisson arrival process versus the number accepted 
Poisson Disk samples is also shown. 



work has demonstrated specialized data structures 
[DH06J for exactly these purposes. In our reference 
implementation, we use a constructive planar geom- 
etry library and approximate disks with polygons for 
simplicity but without a loss of generality. We show 
the performance of our algorithm in terms of samples 
per second, as well as number of samples generated by 
the uniform Poisson arrival process versus accepted 
Poisson Disk samples (see figure [2]). Since the size of 
grid squares is determined by the radius of the PD 
samples, the geometric complexity of the free space 
is 0(1). 

4 Discussion 

We have introduced an algorithm for generating 
Poisson-disk patterns in provable 0(1) time per gen- 
erated sample. Our main insight, compared to re- 
cent O(logiV) per point algorithms, is that rather 
than choosing the location for the next point based 
on area-weighted binary search, we can use an 
area-parameterized exponential distribution to order 



points in time under a uniform Poisson arrival pro- 
cess. While previous algorithms generate each point 
in sequence, with an implicit time linked to their se- 
quential generation, we create many points with ex- 
plicit arrival times and order them (in a local fashion) 
to find those that should be accepted. 
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